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Abstract. Recently, the author and Zia (2010) reported on exact results for 
a far-from-equilibrium system in which two coupled semi-infinite Ising chains at 
temperatures Th and T c , with Th > T c , establish a flux of energy across their junction. 
This paper provides a complete derivation of those results, more explicit expressions for 
the energy flux, and a more detailed characterization of the system at arbitrary T c and 
Th- We consider the two-point correlation functions and the energy flux F(x) between 
each spin, located at integer position x, and its associated heat bath. In the Th — > oo 
limit, the flux F(x) decays exponentially into the cold bath (spins with x = 1, 2, . . .) 
for all T c > and transitions into a power law decay as T c — > 0. We find an asymptotic 
expansion for large x in terms of modified Bessel functions that captures both of these 
behaviors. We perform Monte Carlo simulations that give excellent agreement with 
both the exact and asymptotic results for F{x). The simulations are also used to study 
the system at arbitrary Th and T c . 
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1. Introduction 

The study of non-equilibrium systems has countless applications in many areas, 
including physics, biology, chemistry, and economics. As discussed in a recent review 
of non-equilibrium statistical mechanics pQ, the formulation of a universal description 
of non-equilibrium systems, analogous to the Gibbs framework for equilibium systems, 
has been widely acknowledged as an important goal. Simple models, like the kinetic 
Ising model, provide us with important tools to build intuition and to provide test 
cases for more general treatments. Studies of such simple systems, reviewed in |TJ [2J [3] , 
have greatly contributed to our understanding of non-equilibrium phenomena. One- 
dimensional models are particularly interesting as they are often amenable to analytic 
methods. Moreover, low-dimensional models can be easily simulated on a computer, 
providing us with another tool for exploring non-equilibrium systems. Finally, kinetic 
Ising models are particularly useful as they can be mapped to models of other dynamics 
(e.g., of particles and surfaces) and can therefore characterize a broad class of non- 
equilibrium phenomena with possible experimental realizations [3]. 

We will be interested in Ising models driven out of equilibrium via couplings to heat 
baths with different temperatures. These heat baths set up temperature gradients which 
induce macroscopic energy fluxes in the system. Given the simplicity of the Ising model, 
it is sometimes possible to choose couplings such that analytic results for the energy 
fluxes jl], steady-state correlation functions, corrections to the Boltzmann distribution 
[5j|6], and even the full time-dependent behavior [7] of the system are available. In most 
previous studies of this kind, the systems have a translational symmetry (see [2] for a 
review). For example, a well-studied Ising chain model has alternating spins coupled 
to two heat baths with different temperatures. Infinite range models with multiple 
temperatures have also been considered where all of the spins interact [8] . Finally, there 
have been studies (see [HJOIllEEI]) of quantum spin chains where an energy flux is induced 
directly by an applied field. 

An arguably more realistic way to drive an Ising spin chain out of equilibrium is by 
linking together the ends of two sub-chains, each held at a different temperature. This 
breaks the translational invariance of the system and leads to nontrivial spatial profiles 
for quantities such as the energy flux. Such a localized jump in temperature is found in 
many systems, such as at the interface between the air and a space heater. In the context 
of kinetic Ising models, a possible experimental realization was suggested by Schmuser 
and Schmittmann [6]: Nuclei in a crystalline solid can be prepared at a particular spin 
temperature. Two adjacent domains of nuclei at different spin temperatures might be a 
way to realize the kind of system considered here. Also, some exact analytic results are 
already available for such systems in one dimension (see |T2l [T3] and discussion in jH]) 
where one spin injects energy into an Ising chain via random flipping. In two dimensions, 
models using Kawasaki dynamics (where neighboring spin states are exchanged) reveal 
more interesting features of such driven systems, such as convection cells |15j . In this 
paper, we find exact results for a one-dimensional model where a localized temperature 
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gradient establishes energy fluxes through the system with interesting spatial properties. 

In a recent paper [Hj, the author and Zia presented exact results for an infinite 
kinetic Ising chain with its left half coupled to a hot heat bath and the its right half 
coupled to a colder bath. A sharp temperature gradient is established at the junction 
between the two halves, driving the system far from equilibrium. In this paper, we will 
give a complete account of the brief analysis presented in [H]. We calculate the exact 
expressions for the steady-state two-spin correlation functions and use them to compute 
the flux of energy F(x) between a spin at integer location x and its associated heat 
bath. Intuitively, F{x) describes how energy flows from the hot to cold bath due to the 
temperature gradient. We discuss both the case where the two baths are at arbitrary 
temperatures and the limit in which the hot bath approaches an infinite temperature. 
In the latter limit, we find excellent approximations to the asymptotic behavior of F(x) 
for spins in the cold bath that are far from the junction. 

In section [2] we establish our notation and the microscopic details of the model. 
In section |3] we express the energy fluxes F(x) in terms of the two-point correlation 
functions, which we calculate for arbitrary temperatures for the hot and cold baths. As 
these expressions are rather unwieldy, in section H] we calculate F(x) exactly for the case 
where the hot bath temperatures goes to infinity. For large x, we find that F(x) decays 
exponentially with x for cold bath temperatures T c > 0. As T c — > 0, F(x) decays as a 
power law, with F(x) ~ x~ 3 . We find excellent approximations for F(x) capturing both 
of these large x behaviors in section |5j All of the exact and approximate results for the 
infinitely hot heat bath limit are affirmed by Monte Carlo simulations in section |6j We 
also use the simulations to explore the behavior of F(x) for arbitrary hot and cold bath 
temperatures. We present possible outlooks for future studies and make concluding 
remarks in section [7J 

2. The Model 

We consider a kinetic Ising chain in one dimension with 2N spins. For consistency, 
our notation will be very similar to the notation used in the previous paper [bX]. 
The variables a x = ±1 will denote the two possible values of the spin at site x = 
0, ±1, ±2, . . . , ±N. Although there are many ways to interpret the two values, here we 
will use the language of magnetic spins so that the values ±1 will denote a spin pointing 
up or down, respectively. We begin our analysis by defining a Hamiltonian "H({<t x }) 
for each configuration of Ising spins. For the case of a chain with nearest-neighbor 
interactions with a ferromagnetic constant coupling constant J > 0, the Hamiltonian is 



where we sum over all nearest neighbor pairs (x, y) in the chain. For a chain at 
equilibrium with a single heat bath at temperature T, we can compute a canonical 



N-l 




(1) 



x=-N 



Steady- state properties of coupled Ising chains 



4 



distribution function over spin configurations: P eq ({a x }) = Z~ 1 e~^ n ^ ax ^ , where 
/3 = (A^T) -1 , and Z is the normalization factor or partition function. 

We will be interested in driving our system out of equilibrium by coupling the 
spins to heat baths of different temperatures. Unlike the equilibrium case, we have to 
specify the particular coupling between the heat baths and the spins. The standard 
way to do this is to have the spins transition from one configuration {o~ x } to another 
{a' x } with probability rates ^({a^} — >■ {o~' x }). In the equilibrium case, where all the 
spins are coupled to one temperature, these rates are chosen to satisfy the detailed 
balance condition W({a x } -> K})P eq (R}, t) = W({a' x } -> {a,})P cq ({<}, t). In 
the spirit of the corresponding equilibrium studies, we choose the simple Glauber spin- 
flip dynamics, originally formulated for the equilibrium case [16]. In these dynamics, 
the two configurations {c x } and {o~' x } differ by a single spin flip. The rates W for 
flipping the spins will depend on the temperature of the bath to which the chain is 
coupled. We will generalize these Glauber rates to include a location-dependent heat 
bath temperature T(x) (with its inverse f3(x)). Then, we consider the probability rates 
w x (a x ) = w x (a x — > —cr x ) of flipping a spin at location x (with — N < x < N) given by 

, . 1 



w x [a. 



2At 



0- x \(J x -\ + cr x 



(2) 



2 

where 7(2) = tanh[2/3(a;) J] and (3(x) = [A;bT(x)] _1 , where T(x) is the temperature of 
the heat bath coupled to spin x. The time step At sets the time scale at which the spin 
flipping occurs. For simplicity, we will choose this scale so that At = 1. Since J > 
(so 7 > 0), if the spin at x is anti-aligned with its neighbors, then the flipping rate is 
proportional to 1 + 7(2;) > 1 and to 1 — 7(2) < 1 if it is aligned. Thus, as expected for 
a ferromagnet, the spins tend to align as 7 increases (i.e., temperature decreases). 

We now have to specify what to do at the chain boundaries at x — ±N. One 
choice is to set o~n = o~-n and employ periodic boundary conditions with the same 
flipping rate as in equation Another choice is to allow for open boundary conditions 
in which the boundary spin flipping rates are chosen to satisfy the detailed balance 
condition in the equilibrium case, yielding w±n(o~±n) = [1 — u(±N)a±Na±iy T i]/(2At), 
where u(x) = tanh[/3(x) J] [17]. In the following, we will mostly ignore these boundary 
conditions as we are interested in the behavior of the system as iV — > 00 where we assume 
the effects of the boundaries will be negligible. In a completely rigourous treatment, 
this assumption needs to be checked. In this paper we will check this by comparing our 
analytic results with simulations, which will use the open boundary conditions described 
above. Finally, we will set J = 1 so that all of our energies will be given in units of J. 

When (3(x) is not a constant, the spins can no longer achieve thermal equilibrium 
as the heat baths at each spin x will compete with each other to create a temperature 
gradient, establishing thermal energy fluxes in the system which persist in the steady- 
state. In this case, the detailed balance condition is broken and our steady state 
distribution (which we assume the system goes into as t — > 00) is governed by the 
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time-independent probability distribution 

P*(R}) = P(R},t^oo), (3) 

which is different from the equilibrium distribution P cq ({cr x }). The standard technique 
to deal with this kind of system is to analyze the master equation for the probability 
distribution P({a x },t) of observing a configuration {o~ x } of spins at time t. For the 
non-equilibrium Glauber dynamics we consider, it is given by 



d t P({a x },t) 



1 



N 



w„ 



-a q )P({a x } q , t) - w q (a q )P({a x }, t)} , (4) 



-N 



where the total number of spins is N tot = 2iV + 1 and {o~ x } q is identical to {o~ x }, with 
the exception of a single flip of the q-th spin. The two terms in the summation in 
equation (@J represent the change in P({a x },t) due to transitions into and out of the 
configuration {cr x } via spin flipping at spin q at rate w q (a q ). In this paper, we will be 
interested in the steady-state solution P*({c x }) = P({a x },t — >■ oo) of equation for 
which the left-hand side is equal to zero. 

To fully specify our model, we now consider a localized "junction" at spin x = 
between two Ising chains at different temperatures, as illustrated in figure [U The left 
chain's heat bath will be set to a temperature T(x) = for x = —N, —N + 1, . . . , 0, 
while the right bath will be set to T{x) — T c for x = 1, 2, . . . , N, such that Th > T c . We 
define corresponding j(x) parameters which will satisfy 

for - \ < .;■ < 



(5) 



7 C for < x < N 

and 7 C > 7^. Intuitively, we expect the temperature gradient to induce a flow of heat 




t 



Ifcu 



-1 -3-2-1 12 3 4 



t 



5 



x 



Figure 1. A schematic of two adjacent kinetic Ising spin chains coupled to heat baths 
at different temperatures Th and T c , with Th > T c . The red arrows show the flow of 
heat F(x) between the heat baths and the chains. Heat flows across the junction from 
the hot bath and into the cold one. As \x\ increases, the heat flow will decrease as the 
spins approach equilibrium with their respective heat baths away from the junction. 
This is illustrated by the decreasing size of the red arrows. We expect that the heat 
flux decays faster in the hot bath, as discussed in the main text (see section [6]). 



from the hot bath into the cold bath. The most dramatic non-equilibrium behavior 
will occur at the junction between the chains at location x — 0. Far away from this 
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junction, i.e. for x — > ±00 in the thermodynamic limit (N — > 00), we expect that 
the chains are at equilibrium with their respective baths and the heat flow decays to 
zero. An illustration of the model and the expected qualitative behavior are given in 
figure [TJ Our objective now is to compute how this energy flux F(x) decays away from 
the junction as a function of the spin location x in the stationary state characterized by 

p.(W). 

3. Two Heat Baths at Arbitrary Temperatures 

In order to find the net energy flux F(x) at each spin, we recognize that all contributions 
to F(x) will come from spin flips at spin x. Then, given our Hamiltonian in ([1]), we see 
that the energy change AE(x) gained by the heat bath (or lost by the chain) due to a 
single spin flip at location x is given in units of J as 

AE(x) = H({a z }) - H({a z } x ) = -2a x {a x . x + a x+1 ), (6) 

where a x is the value of the spin at x before the flip. We also know the flip rate u x (o~ x ) 
at spin x. Therefore, the net heat flux F(x) must be given by [I] 

F{x) = (u x (a x )AE(x)) 

= 7(x)(l + (o- x -io- x+1 )) - (a x a x+ i) - (cr x o- x -i) . (7) 

As we are interested in the average solution in the stationary state, the bracket averages 
(. . .) are with respect to the stationary distribution P*({cr x }). We can think about F(x) 
as a flow of heat from the spin chain into the heat bath. We have F(x) < for heat 
flowing into the chain and F(x) > for flow out of the chain. 

We see that to compute F(x), we have to consider the two point correlation 
functions 

(o- x O-y) = y^ j Q- x Q'yP*{{0'z}), (8) 

{».} 

for which we define a convenient notation (x, y) = (o~ x o~ y ). Using a calculation 
completely analogous to the one done by Glauber for the equilibrium case [IS] , it is 
possible to derive a difference equation for (x, y) from our master equation in (jl]). First, 
we set the left-hand side of equation (Tj0) to zero to get an equation for P*({cr z }). Then 
we multiply this equation by a x a y and sum over all configurations {cr 2 }. This gives us an 
equation for the two-point correlation functions (x, y), which, after some manipulations, 
we can write as 

7 (x) ((x + l,y) + (x- 1, y)) + 7(1/) ((x, y + 1) + (x, y - 1)) - 4 (x, y) = 

[ 7 (z)«^ + 7 (y)aj + 2 ( 7 (x) + 7 (y) - 2)] (x, y) = 0, (9) 

where 5 X is a second order difference operator which acts on any function f(x) as 
Slf(x) = f(x + 1) + f(x — 1) — 2f(x). Notice that equation ([91) has the form of 
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an anisotropic discrete Helmholtz equation, which can be solved using Green's function 
techniques. To do this, we require appropriate boundary conditions. Given the definition 
of the two-point correlation function in (JSJ, we must have the two boundary conditions 
(BCs) 

(x, x) = 1 

V ' ' . (10) 

(x,y) = (y,x) 

Another boundary condition comes from the behavior of the system for large separations 
\x — y\. Namely, we expect that as \x — y\ — > oo, the spins become uncorrelated so that 
(x,y) — > (o~ x ) (o-y) = 0, since there is no spontaneous magnetization in one dimension 
and (a x ) = for any choice of < j(x) < 1 and x, as discussed in [TJ]. Since we are 
dealing with a finite system for now, we will have to implement this boundary condition 
by enforcing it exactly at the boundary spins at x = ±iV. So, we will also have the BCs 

(x,±N)=0 (11) 

for all x ^ ±N. 

The second condition in ( |T0|) implies that we can now just consider x > y for the 
purposes of calculating the correlation functions. Then, given our definition of 7 (x) in 
([5]), we can identify three regions in the (x, y) plane on which we must solve equation 
([9]) with the appropriate BCs. These regions are illustrated in figure El Substituting 
the appropriate values for j(x) from ([5]) into equation (J9j), we find the equations 

[lh$l + Ih^y + 4 (7^ - 1)] (x, y) = in R h : y < 0, x < 

M + lh% + 2 (7c + lh - 2)] (x, y) = in R hc : y < 0, x > , (12) 
[lc5 2 x + l c 5 2 y + 4 ( 7c - 1)] (x, y) = in R c : y > 0, x > 

the solutions to which must be matched along the red lines x = and y = shown in 
figure El 

To find the solutions to equation ([9]) in the various regions, we look at the Green's 
function G(x,y; £,77) in each region Rh,hc,c that satisfies 

V x>y G(x, y,tr])= [l{x)5 2 x + 7 (y)<jJ + 2 ( 7 (x) + 7 (y) - 2)] G(x, y;Z, V ) 

= 5 x ^5y tV (13) 

where 7(x) is given in (|5]), 5 x , y are Kronecker delta functions, and T> XiV is a convenient 
notation for the discrete operator in the square brackets. Our Green's function is chosen 
to satisfy G(x, y ; £, rf) = along the red, green, and blue lines in figure [13l Performing 
a Fourier transform, we find a solution for G(x, y; £,■)]) for values of (x, y) and (£, rj) in 
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Figure 2. The regions Rhc, and i? c represent correlations between spins coupled 
to the hot, hot and cold, and cold bath, respectively. We want to solve the equations 
in (Tl2j) for these correlations. The dashed green lines (y — —N and x — N) correspond 
to the BCs given in (fTT|) . We must enforce the first condition in (fT0|) on the dotted 
blue line {y — x). Finally, we have to ensure our solutions match along the boundaries 
between the regions indicated by the solid red lines (x = and y = 0). 



the three regions of interest: 



E E 

m=-N n=—N 
N 



U mjn (x,y)U m<n (£,r)) 



EE- 

m=0 n=— 
N m 

EE 



4 - 2^ h [cos (irm/N) + cos(irn/N)} 



in R h 



4 — 2 / y c cos(irm/N) — 2^^ cos(7rn/iV) 
U m , n {x,y)U m!ri (^,r]) 



in i? 



/if 



(14) 



in 



,^1^ 4 - 2 7c [cos(7rm/iV) + cos(v™/iV)] 

where we take the linear combinations of Fourier eigenfunctions which vanish on the 
region boundaries and are normalized over the regions: 



U m ,n y) 



2 






sin 


N 




2 




— sin 


TV 





firmx 



sin 



sin 



( 



nny 
nny 



— sin 



rrmy 
~N~ 



sin 

\ N J\ 



(15) 



Using this Green's function, we are now able to solve for the two-point correlation 
functions. 
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We find the correlation function (x, y) by exploiting a discrete Green's theorem. 
Namely, if we multiply equation f lT3|) by (£, 77) and sum over all values of £ and 77 in the 
region we are interested in, we find that 

(x, y)= Y, ^ ;Z,V) - G(x, y;Z,ri) V itV (£, 77}] , (16) 

where we recognize that T>^ n (£, 77) = 0. One can show that the summations over all 
(£, 77) G Rh,hc,c m (fTBI) reduce to summations just over the boundaries dRh t hc,c of each 
region. In particular, for the three regions we are interested in, we have the solutions 
(using the more compact notation G x>y (£, rj) = G(x, y ; £, 77)) 




rj=-N 

N 

ri=-N g=l 
N 



in Rh 



in i? 



(17) 



in R r 



where we have simplified some boundary terms by applying the boundary conditions 
(Eqs. [TOl [11]). Finally, we must find the functions (£,0) and (0,77) that match the 
solutions in the three different regions and are consistent with Eq [9l We now have a 
complete description of the general solution, but finding an analytical form for (x, y) is 
cumbersome. To make progress, we will now look at the energy flux through the chain 
in the case that the hot bath is infinitely hot. This limit will allow us to get a much 
more detailed understanding of the behavior of this model. 



4. Hot Heat Bath at Infinite Temperature 

Let us now set 7 = j c , T = T c , etc. We then set 7^ = to keep the hot bath at infinite 
temperature. We immediately see from ( [12]) that all of the correlations in region Rh 
vanish. In region Rh c all the correlations vanish, except for the ones at the boundary: 
(x, 0) for x > 0. This is intuitive because, for Th — >■ 00, all of the spins in the hot bath 
flip randomly. Thus, any spin in that bath is equally likely to point up or down and 
cannot be correlated with a non-neighboring spin in the cold bath. Explicitly, we see 
from ( fl2i) that for 7^ = 0, the correlation function (x,y) satisfies 

[76I + 2 (7 - 2)] (x, y) = 7 ((x - 1, y) + (x + 1, y)) - 4(x, y) = (18) 

for all y < 0. This is a homogeneous, linear second order difference equation which we 
can easily solve using standard techniques. Namely, we try to find the solution of the 
form (x, y) = Ar x . Substituting this ansatz into equation ()18p and solving for r, we find 
that there are two possible values of r so that the general solution for (x, y) is given by 

(x,y) = A 1 u x + A 2 u- x , (19) 



Steady- state properties of coupled Ising chains 



10 



where A\2 are constants set by the BCs and 



2 / 4 

u=--J--l. (20) 
7 V 7 

Then, since all of the correlations in Rh vanish, we know that (0, y) = for all 
y = —N + 1, . . . , —1. This means that we must have A\ = —A2 = B /2 so that our 
general result in Rh c becomes (x, y) = B sinh(x In a)). We now let N — > 00 and require 
that that (x, y) — > as x — > 00 for any y. The only way this can happen is if B — 0. 
This means that the (x, y) vanishes in region Rh c , as well. The energy fluxes for x < 
must also vanish since we will have, via ([7]), 

F(x) = -[(x,x + l) + (x,x-l)] = for all x < 0. (21) 

We now compute the correlations in region R c and the boundary function (x, 0) for 
x > 1. The latter function satisfies equation ( fl~8l) with y = 0. Therefore, we have the 
same general solution given in (Tl9|) . This time, our boundary conditions are (0,0) = 1 
and (x, 0) — > as x — >■ 00. Since < a) < 1, these boundary conditions set A 2 = and 
A\ — 1 in ([19]) , yielding the solution 

(ar, 0) = w a for x > 0. (22) 

Consequently, at the x = spin, we have an energy flux equal to 

F(x = 0) = -u. (23) 

We have now shown that all of the flux from the hot bath into the spin chain is 
located at the boundary spin at x = 0. The negative sign in the expression for F(x = 0) 
means that, as expected, the heat flows into the chain from the infinitely hot bath. 
With equation ( TT71) in region R c , we now have a complete solution to our correlation 
functions (x,y) and, consequently, F(x) for all x > 0. However, the resulting expression 
for F(x) is cumbersome. To find a simpler solution, we now recognize that the Helmholtz 
equation in R c (see ([121 ) has the same form as for an equilibrium Ising chain (a chain 
of spins all coupled to a single heat bath with a parameter 7). The only difference 
between the equilibrium solution (x,y) eq and the non-equilibrium one (x, y) must be in 
the BCs. The (x, x) = 1 BC is the same in both cases, so the difference must come 
from the (x, 0) correlation we computed in ( 122]) . In the equilibrium case, we can find 
this function by solving equation ([9]), setting j(x) = 7(2/) = 7. We find the standard 
exponentially decaying result [H] 



(x, y ) eq = Jy-^ 



.7 V 7 2 1 



\y-x\ 



\y-x\/€e 



(24) 



where we have the equilibrium correlation length £ eq = — (lnw) -1 and we have used 
the relation uj = 7 _1 — a/7~ 2 — 1 for oj = tanh/3 and 7 = tanh(2/3). Notice that the 
equilibrium and non-equilibrium solutions for (£, 0) are related by a replacement of 7 
with 7/2, the average of the parameters of the hot and cold chain: and 7. 

Given the close relationship between the non-equilibrium and equilibrium cases, it 
is convenient to study the deviation A(x,y) = (x,y) — (x,y) cq . Both (x,y) and (x,y) eq 
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satisfy ( I12p in R c with 7 = 7 C and the linearity of the Helmholtz equation allows us to 
use (j!7p to solve for A(x, y). The only difference is in the boundary condition along the 
y = line. There we must have 

A(x, 0) = u x - u x . (25) 

Thus, the solution to A(x,y) in R c is 



N 



A(x, y) = 7 ^2 y ; C ? !) - ^] for < y < x. 

5=0 

Substituting the expression for G(x, y ; £, 1) in equation ( fl4l) into (l26l) . we find that 

TV m-1 TT I \ 

A( 1= 7yy UmA^y) 

lX ' W AT 2^ 2 - 7 [cos(vrm/iV) + cos(vrn/iV)] 



(26) 



E(- £ -- £ ) 

.5=0 



sm 



sm 



/7rn\ /7rm\ 



\ N J 



sin 



irnC, 



N m-1 rr 

TV 2^ 2^ 



(2, y) sin(7rm/iV) sin(7rn/iV) 



m=l ?i= 



J 2 — 7[cos(7rm/A^) + cos(7m/iV)] 



1 — 2a; cosf 



A? 



7,2 



+ 



l-2wcos(^) +w 2 l-2wcos( T ^) +w 2 l-2wcos(^) 



(27) 



where the sum over £ was performed by expanding the sin(7rm£/iV), sin(7rn£/iV) terms 
in exponentials and summing the resultant geometric series. We have also made the 
approximation that u N ~ u N ~ for large N. This approximation will be exact in the 
iV — > 00 limit for any < 7 < 1. We now move to the calculation of the flux F(x) for 
x > 0. 

Since F(x) = at equilibrium, we can replace (. . .) with (. . .) eq in (J7J) to get 
a combination of equilibrium correlation functions that must equal zero. Then, we 
subtract this combination from the right-hand side of ([7]) to find 

F(x) = 7A(x + 1, x - 1) - A(x + 1, x) - A(x, x-l). (28) 

Since ( 128]) includes the functions A(x, x — 1) and A(x+1, x) and we impose the boundary 
condition in (l25l) . we will have to treat the x = 1 case separately. For now, suppose 
x > 1. Then, substituting ( 12"TI) and (ITS]) into equation (T2"8"|) . and using the identities 
uj = 27 _1 — a/47~ 2 — 1 and u; = 7 _1 — a/7~ 2 — 1 , we get 

2 N m 



F(x) 



EE 



m=l n=l 



sin (™*) cos sin (^) (l - 7 cos ($)) 



X 



2 -7 [cos + COS (f)] 

3- 7 cos(f) - 7 cos(^) 



— n -H- m 



[2 -7 cos (f)] [I-7 cos (sp)] 



— n -H- m 



sm 



sm 



\ iV / \ iV / ' 



(29) 
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where the n -H- m term is just the first term in each bracket with n and m exchanged. 
Moving to the N — > oo limit, we define the continuous variables 



k 



Tim 



and 



7rn 



(30) 



The summations in f )29|) become integrals as we can identify them as regular Riemann 
sums with step size tt/N. Also, notice that the summand in (12"§|) is symmetric under 
the exchange of m and n, vanishes for m = n, and is even in both m and n. So, moving 
to the N — > oo limit, making use of the symmetry arguments, and simplifying some of 
the terms gives us 

ship (1 — 7 cos k) — sin k (1 — 7 cos p) 



F(x) = Re 



V 
Ain 2 



dkdpe 



i(p+k)x 



— TV J —TV 



X 



1 



1 



1 



2 — 7(cos k + cosp) 
1 



sin k sin p 



(31) 



2 — 7 cos k 2 — 7 cos p 1 — 7 cos p 1 — 7 cos k 
where we have written our expression as the real part of a complex function that will 
aid us in what follows. 

To make progress calculating the integral in (I3"T|) . it is convenient to change variables 
to s = k+p and q = k—p. Since our region of integration is k G (— 7T, 71") and p G (— 7T, 7r), 
we can let g range from (— 2tt, 2n) and have s G (|g| — 27r, 27r — |g|). Changing to these 
variables and applying some trigonometric identities yields 



F(x) = Re 



2tt 



27T-IJ 



H x (s,q) ds dq 



-2tt 



(32) 



where 



H x (s,q) 



[cos 2 (f) - cos 2 (I)] [7 cos (I) — COS (|)] 
l- 7 cos(|) cos(f) 

(-l) a sin(f)sin 2 (f) 



E 



(33) 

^— ^ 2a 2 — 47a cos (|) cos (I) + 7 2 (cos s + cos g) 

We have also recognized that the real part of the integrand in ( 132]) is even in q so that 
we can change our domain of integration to q G (0, 2ir) and s G (2ir — q,q — 2n). We first 
compute the integral over s in ( 1321) via a rectangular contour integration in the complex 
s plane. The contour C is shown in figure |3j Cauchy's residue theorem tells us that our 
integral of interest (the one over Co in figure |3]) is given by 

/ ds H x (s, q) = 2m V Res - / H x (s, q) ds 

Jq-2-K - x JCi 



H x (s,q) ds 



c 2 



H x (s,q)ds, (34) 



C*3 



where the first term represents the contribution from the poles inside C located at s* , 
with residues Res[sl' 7 ' ) ] . We shall see in the following that there will be five poles. 
Before we move on to these poles, consider the integrals over the various parts of the 
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contour C. First, the contribution from C2 vanishes since the factor e lsx = e ltx ~ TX will 
go exponentially to zero as r — > 00. The integral over C\ is a bit trickier to analyze and 
is given by 

p /»oo /*oo 

/ dsH x (s,q) = i dtH x (27r-q + it,q)= dth Cl (t,q). (35) 
Jd Jo Jo 

Using ( l33l) . we can show that the real part of hc\(q,t) satisfies, for integers x > 0, 

Re [h Cl (t, 7T + q')] = -Re [h Cl (t, vr - q')] , (36) 

for all q' G (— 7r, 7r) and t > 0. Thus, the contribution from C\ to F{x) must vanish 
since we will have 

/ / H(s,q)dsdq \ = / / Re [h Cl (t, vr + q')] dq' dt = 0. (37) 

The contribution from C3 vanishes by an analogous argument. All that remains is to 
evaluate the contributions from the poles. 



Im s * 

T — > OO 



c 2 



Czl 



1 s? ] n ln[«- 3 ]" 



si 4 ^ X lu [ w 2 ]" 



— h- 

q-2n -q 



s {1) 



(3) 



(5) 



Re s 
♦ 



P Co g 2?r-g 



Figure 3. Here we denote the contour C = Co + Ci + C2 + C3 in the complex plane 
that we require to perform the integration over s in (|32]). We will let r — > 00 so that 
the contour becomes a semi-infinite strip. The red crosses denote the locations of the 
poles that are picked up by the contour. The locations of the poles si^ (i = 1, . . . , 5) 
picked up by C. 



Notice that the functions in the numerators and denominators in the expression 
for H(s, q) in fl33|) are all analytic over the entire complex s plane. Thus, the poles will 
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come from the zeros of the denominators, i.e., the solutions to the equations 
l- 7 cos(£)cos(f)=0 

2a 2 — 47a cos cos + 7 2 (cos s + cos q) — a = 1, 2 
There are many possible such solutions, but only five will be picked up by the contour 



(38) 



C, as shown in figure [3) These five poles are located at s* 



('■) 



,5), given by 



sjf' = % 2 In 
>f ' 3) = ±g-i21nw 



- lsec (|) + J I' 2 sec2 (| 



(39) 



v 4 4 ' 5) 



±q — i 2 In Co 

These poles are only picked up by the contour C when < q < n: the real parts of s* 
for i — 2, 3, 4, 5 will not be between q — 2tt and 2n — q when q > ir, and s* is pushed 
off to +200 as g — > 7T. So, we now can change our limits of integration for the integral 
over q: 



F(x) = - 



5 pjY 

2vr V / Im [Res (s[ j) )] dq. 
3=1 Jo 



(40) 



We now compute the residues at each pole. We deal with the i = 2, 3 poles first. The 
residues are, after much simplification, 

L J An 2 (co + ur) 

with the upper (lower) sign corresponding to the residue at s^ (s^). Picking up the 
imaginary part and performing the integrations yields the same result for both poles (as 
we expect since Re[G(s,q)} is odd in the s variable): 

,2*-i ( w 2 _ 1) [1+ aj + w 2(i _ x )] sinjvrx] 



Im{Res [sl 2 ' 3) ]} dg 



4tt 2 (x 2 - 1) (1 + u 2 



uj(1-uj 2 ) 7 (1 -a; 2 ) 

Oil — o Gal > 



(42) 



4tt(1 + u; 2 )~^ 8vr 
where we see that when x is an integer, this contribution vanishes for all integer x > 1 
and has a finite contribution at x = 1. Since we are dealing with the x > 1 case, we can 
ignore this contribution. It will, however, turn out to be useful later. 

We find a similar situation for the i = 4, 5 poles, with slightly more complicated 
integrals over q. Here, the residues are 

r (4 5)1 _ ie^+V (eg - 1) u 2x {u 2 - e^ 2iq ) [e^(l - 3u 2 ) + u 2 {u 2 - 3)] 

[S * ' J ~ 4tt 2 (e^ + u 2 )u(l + u 2 ) ' { } 

where the upper (lower) sign refers to the residue at si^ (s* ). To perform the integral 
of this contribution of q, we first expand f !43p out in terms proportional to e iqc for some 
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we find that all of the terms involve integrals 



constant c. For the integral of Res[s 
of the form 

4(C) 



dg 



e -2iq + u 



~,2 



iuj 2 



dz 



Z 2 + LO 2 



(44) 



where C > —2 is a constant and we have made the variable substitution z = e* 9 so that 
our new path of integration is S + , the unit semi-circle in the top half of the complex 
z plane, traversed clockwise. There are no poles between S + and the real axis because 
the poles are located at z* = ±iu~ 2 and we know that uj~ 2 > 1 for all < 7 < 1. Thus, 
deforming the contour over S + onto the real axis yields 



4(0 



1 z (+l 

l 1 + bJ l Z l 



du 



u 



C/2 



1 + LU 2 U 



= <[l + (-l)«+ 1 ]F (C), 
where we have made the substitution u = z 2 , defined 

1 



mo = 



C + 2 



c c 



-u 2 



(45) 



(46) 



and recognized in ( 145]) a standard integral representation for the hypergeometric function 
F(a, /3,7;z) (see e.g., equation 3.197 3 in [IS]). 
Using (145]) to integrate Res[sl 4 ' ) ], we find that 



dgRes [ S < 4 >] 



~,2x ( A-kx 



+ 1) 



4tt 2 w(1 +cu 2 ) 



4tT 2 Cu(1+£ 2 



(1 - 3w 2 ) F (a; - 3) + 2w 2 (w 2 + 1) F c 
+w 4 (w 2 -3)^(x+l) 



1) 



(47) 



The integral we encounter for the Resfsl 5 " 1 ] term is very similar. However, we do not 
have to calculate anything for the residue at sf^ because we can recognize from ( 143]) 
that the imaginary parts of the residues at si 4 '* and si 5 ^ are the same. Thus, we just have 
to examine ( l47j) to find the contribution from both residues. We find that for integer 
x > 1, these integrals give us purely real contributions. The only imaginary contribution 
will come from the limit x — > 1, where the first term on the right-hand side of equation 
( ]47j) picks up a factor of — in since (e i7rx + l)/(x — 1) — » —iir when x — > 1. Therefore, 



Im{Res [4 4,5) ]} dg 



£J (3u 2 - 1) 



7(3w 2 - 1) 



5xi- 



(4£ 



4tt(1 + w 2 ) 16vr 

We now conclude that the only contribution to F(x) for x > 1 must come from 
the sl 1 ^ residue. The residue turns out to be pure imaginary and substituting it into 
equation ( 140]) gives us 



F(x) 



4 

^7 Jo 



(7*7) 1 - V(jnr 



1 



2a- 



1 



7V) 



[i + ^2 {v 2 _ !)] yy-r^ 



dr/ for a; > 1, (49) 
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where we changed variables to 77 = cos(g/2). Let us finally deal with the x = 1 value of 
the flux. We substitute (1251) into equation (1281) to find that 

F(x = l) = 7(w 2 -w 2 )-w + a;-A(2,l). (50) 

Notice that the x — > 1 limit of the expression for F(x) in (1291) coincides with the equation 
for — A(2, 1) in ( 1271) . So, we are now able to use the contributions we found in equations 
(132]) and (iE]) to calculate A(2, 1). We find 

(iv)' 1 



Fix 



7(o, 2 



0J ) — u + u + 



7V) 



+ 



7 (1 — 3cj 



~,2\ 



^7 JO 

7(1 



[1+ 7 V (t? 2 - 1)] v 7 ! 3 ? 



d?7 



4 

ti-7 Jo 



(777) 1 - \firn) 



-2 



7V) 



dry. 



(51) 



[I + 7V (?7 2 - 1)] a/1 - ?7 2 

We now see that the behavior of F(x) for all x > 1 is governed by the residue at s^ . 
This was previously argued for just x > 1 in [13] . but we have now carefully calculated 
F(x) for all x > 0. In summary, we have the exact expression 



F(x) 



4(1 - 5 x0 ) 

7T7 



(777) 1 - ^(777) - 



1 



2.r 



7V) 



d?7 - ( 52 ) 



[1 + 7 2 77 2 (772 - 1)] a/1 - 77 2 
for all integers x > 0. 

The flux of energy must be conserved globally so that the flux into the system 
from the hot bath (—&) must be compensated by the total flux out of the cold bath 
Ylx>o F{ x )- Since the expression in the brackets [. . .] 2x in fl52|) is in the interval (0, 1) for 
any 7, 77 G (0, 1), we can perform the summation over all x to find the amusing integral 
identity 

/ , x 2 

1 00 o yl (1 



x=l 



7V) 



7T7W7o 0^ [1 + ^2 - 1)] 7^' + _ 1 



1 



(53) 



for all 7 G (0, 1). This identity can be evaluated numerically as a kind of check of our 
exact solution. 



5. Energy Flux Asymptotics 



Looking at the large x behavior of F(x) allows us to characterize the most important 
physical features of the model. We shall see that the asymptotics tell us about a cross- 
over in F(x) from an exponential to a power law decay as T — > 0. To find these large x 
asymptotics of F(x), we have to examine the integral in ( 152]) 



If(x) 



eX P \ ~ 2x ln 



[7 cos 9} 1 + \/ [7 cos 9} — 1 



7 2 cos 4 1 



(7/2) W (20) 



■d0. 



(54) 
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where we have restored the original integration over q and changed variables 9 = q/2. 
We know the main contribution to Ip will come from the minimum of the function 

1 



f(9) = In 



+ 



1 



(55) 



7 cos 9 ' V 7 2 cos2 6 

This minimum occurs at 9 = 0. In order to be able to apply Watson's lemma [20] to 
our integral, we change variables to 

s = 2/(0) + 21nw, (56) 

which will satisfy s G (0, oo) whenever 9 G (0,7r/2). Inverting equation (1561) and 
performing this change of variable, equation ( 1541) becomes 



7 F (x) 



where 
W(s) = 

and 



„2a: In c 



7^) 



dsW(s)e" 



2 [1 + Th£{rfi-l)]y/T=ti 



r 



cosh 



+ sinh 



Vs = 



cosh 



a/1 — 7 2 sinh 



(57) 



(5f 



(59) 



The expression for W(s) is extremely unwieldy, but we now see that the integral in ( 1571) 
is in an appropriate form for the application of Watson's lemma. Thus, we only need 
to worry about the expansion of W(s) around s = 0. At s = 0, we see from ( )58l) that 
W(s) develops a square root singularity s _1//2 due to the a/1 — ^1 term, which vanishes 
like a/1 — r]l = (1 — 7 2 ) 1 / 4 y/s(l + O(s)) as s — > 0. If we extract this singularity, the 
rest of W(s) is analytic in s, which we can expand as a Taylor series around s = with 
coefficients a n . This yields, 



W(s) 



;i-7 



n=0 



2A/i 



1 



(87 4 - 297 2 + 2) s 



+ 0(s' 



(60) 



This expansion and subsequent ones were performed with some assistance from the 
computer algebra system Mathematica 8.0. The expansion in (160!) becomes invalid in 
the 7 — > 1 limit, as the higher order terms a n in (|60l) become larger than the lower order 
ones, becoming divergent for n > 3. We will deal with this problem later. For now, we 
may conclude by Watson's lemma |20j that for 1 — 7 sufficiently small, If(x) admits an 
asymptotic expansion given by 

" a n r(n + l/2) 



If(x) 



2x In a) 



E 

71=0 



x n+l/2 



So, to second order, the flux is 



F(x) 



2(1 



7 2 ) 5/ V 3; 



7a/vtx 



1 + 



29 7 2 _ 8 7 4 - 2 
16xa/1 — 7 2 



+ 



x 



r 



-2 



(61) 



(62) 
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The u 2x = e 2xlnuJ term in (1621) tells us that the energy flux F(x) decays exponentially 



as e x ^, with a decay length equal to 

e= 1 



21nw 



2 In (j- 1 



'1 



-2 



2 ' 



(63) 



which we recognize as half the equilibrium correlation length £ eq (see (|Hj) ). Notice that 
the quantity x\Jl — 7 2 has to be /arge for our expansion in (162]) to remain valid. Since 
this is impossible for 7 — > 1 (T c — )■ 0), this approach is invalid in this limit, and we must 
use another method. 

To see what happens as the cold bath temperature drops to zero, we first take 
the 7—7-1 limit in the expression for W(s) in fl58|) . We then find a different analytic 
structure in the s variable with no square root singularities: 

2(3 + cosh s) sinh (§) tanh (f ) 



W(s) 



1 7— >i 



_ s* sr 
~ T ~ 24 
Watson's lemma then tells us that 



7 + cosh(2s 
179s" 



E<>« 

n=0 



,2n 



23040 



+ C(/ 



. 00 

—E 

n=0 



&n(2n)! 
l+2n 



71X' J 



TCX" 



+ 0(x- 7 ). 



(64) 



(65) 



(66) 



We see that _F(x) transitions into a power-law decay in the 7—^0 limit. This behavior 
is reminiscent of an equilibrium system at a critical point where we find long range 
correlations. The particular leading order behavior, F(x) ~ x~ 3 , in this non-equilibrium 
system is interesting and does not seem to have a simple explanation. 

To capture both the power law decay of F(x) in the 7 — > 1 limit and the exponential 
decay for 7 < 1, we can no longer apply the analysis utilizing Watson's lemma described 
above. Instead, we want to set up an expansion that captures both the xy/ 1 — 7 2 3> 1 
regime, where F(x) exhibits an exponential decay, and the x^/l — j 2 1 regime, where 
F(x) has a power law decay. To do this, we go back to equation (l54"j) and rewrite it in 
the modified form 



If(x) - 
where 

and 

g{y) = 



2x 2xy/l-j 2 



f 



COS 



2/J 



e 2x 9 [cos( g /2)] exp 



-xy/A(l - 7 2 ) + q 2 ] dq, (67) 



f(y) 



2„,4 



1-72/ 



i 2 y 2 (y 2 



a/(1 - 7 2 ) + [acos 



-In 



(iy) 1 + VJjyf 



•A 



7 



lnu. 



(68) 



(69) 



As in the Laplace method, we now argue that the main contribution to the integral in 
(|67|) will come from the region q pa 0. So, we now perform the expansion by analogy 
with the Laplace method by expanding the term f{y)e 2x9 ^ for y = cos(g/2) in a power 
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series in q and extending our range of integration to ( G (0, oo). From the expressions 
in (168]) and ( 169]) . we see that only the even powers of q will contribute to the expansion. 
To facilitate this expansion, we first define the constant ( 2 = 4(1 — 7 2 ) and calculate 
the integral 

/*00 /*0O 

J g 2n e -,V(V d? = C 2n+1 J z ^2 _ ^ ^ 

C 2n d 
/tv dx 



2 n (x()- n T[n + -)K n (xC) 



(2n-l)!! (20" 



7' 



n+l 



AT. 



X' 



n+l 



2xa/T 



7' 



(70) 



where n > is an integer, n!! = n(n — 2) . . . for any positive integer n (with 
0!! = ( — 1)!! = 1), and we have recognized an integral representation of the modified 
Bessel function (see e.g., equation 3.387 6 in [19J). Using flTOj) . we may now set up our 
asymptotic expansion of the integral in ( 1671) of the form described in [14"] : 



where 



If(x] 



B„ 



n=0 

(2n-l)!![4(l-7 2 )] n+1 / 2 d 2n 



f 



cos 



2xg[cos(q/2)] 



(71) 



• (72) 



9=0 



2(2n)! dq 2n 

These coefficients can be computed exactly for any n. We find that if we keep the first 
two orders in 1/x, we have 

Au 2x (l - 7 2 )e 2x V^~ 2 



F{x) 



7T7 

7 2 (3 -7 2 ) 
2x 



r 



K 2 [2xy/l-r 



5(1 -7 2 ) 
16x 



K 3 [2x^/T^)+0(x- 2 ) 



.(73) 



One can easily check that the expansion in (1731) captures both terms in equation 
for (x ^> 1 and reduces, for (x 1, to 2x~ 3 /n. We conclude that this Bessel function 
expansion faithfully captures the behavior of F(x) in both of the regimes of interest. In 
the next section, we will compare these asymptotic results with both simulations and 
the exact expression for F(x) in ( |52l) . 



6. Simulations 

We will now verify our solutions for F(x) via a Monte Carlo routine. In this routine, 
we initialize our Ising chain of 2N + 1 spins at locations —N < x < N in a random 
configuration and then evolve it according to the master equation rates. This is done 
by choosing any spin a q with equal probability and flipping that spin with probability 
Atw q (a q ). We then advance the time t by t — > t + At and repeat the process. We choose 
open boundary conditions at the ends. In the Th — > oo limit, the spins coupled to the 
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hot bath flip with probability 1/2, regardless of the states of their neighbors. Therefore, 
we can perform our simulation in this limit on just the N spins coupled to the cold bath 
and the single infinitely hot spin at x = 0. 
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Figure 4. A log-log plot of time time dependence of both the mean AE(x, t) and the 

2 

variance AE(x, t) 2 — AE(x, t) for 7 = 0.8 and x = 5. The blue lines show a fit of the 
data to a linear function and indicate that both the mean and variance grow linearly 
at large t. The y-intercepts of the lines on the log-log plot can be used to estimate the 
slope a of the fit and, thus, $5 and V5, as in ([74]) . 



To record the flux F(x) at each spin, we have to keep track of the spin flips at 
each spin as each can contribute a change in energy of ±4, depending on whether the 
spin is aligned or anti-aligned with its neighbors. We compute the average change in 
energy AE(x,t) over about 3 x 10 6 independent simulation runs after t = 2 n Monte 
Carlo updates at each spin (for n = 0, 1, 2, . . . , 20). We tested values of 7 between 0.5 
and 1. For the smaller values of 7, the sharp exponential decay of F(x) required more 
simulation runs (up to 1.2 x 10 7 ) to get good statistics at large values of x. We can 
think about these Monte Carlo updates at each spin as inducing a random walk in the 
variable AE(x,t), with "step size" 4. 

In an equilibrium case, we expect that the average AE(x,t) is zero, as we have no 

net flux F(x). However, in our case, the presence of the infinite temperature spin at 

x = induces a bias in the AE(x,t) dynamics. The non-zero value of F(x) induces a 

uniform drift and, thus, AE(x, t) is non-zero and grows linearly in time after the system 

achieves a steady-state at long times t, as shown in the log-log plot in figure HI Notice 

that the data points for AE(x = 5,t) versus t fall perfectly on a line with unity slope 

on a log-log plot for large t, implying linear growth. We also see in the figure that the 

2 

variance, AE 2 (x, t) — AE(x, t) , also increases linearly in t, just as we would expect for 

a random walk. In summary, we use these simulation results to estimate the two slopes 
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V x and <& x , where 



AE(x,t) « $ x £ 



A£(a;, t) 2 - AE(x, t) « ^ 



(74) 



As illustrated by (j7]), the rate of increase of AE(x, t) in the steady-state of the system is 
equal to the energy flux F(x). Thus, our flux F(x) is estimated by F(x) $3.. We can 
also approximate an error in this estimate by utilizing properties of the random walk. 
The standard error of the mean AE(x,t) is given by 



SE [AE(x,t) 



AE(x,t) 2 - AE(x,t) 
N 

1 » runs 



(75) 



where iV runs is the number of runs over which we average. Then, given the definition 
of in ( 1741) . we see from ( J75|) that an estimate of the error ±AF(x) in F(x) ~ Q x is 
given by 



AF{x) 



Vr 



(76) 



N t 

1 ' runs b 

where t will be the number of time steps used to determine our slopes $> x and V x . 

An easy initial check of our simulation is to plot the flux at the zeroth spin F(x = 0) 
versus the cold bath parameter 7. We already argued in (I23p that we have the exact 
solution F{x = 0) = —u = a/47 -2 — 1 — 27 _1 in units of J. We see in figure |5] that the 
simulation results agree perfectly with this exact result for a wide range of 7. Notice 
that as the cold bath temperature T — > 0, we have 7 — > 1 and the flux increases in 
magnitude. This makes sense as we increase the temperature gradient between the 
chains as we decrease T, forcing more energy through the junction between the chains. 




-0.08 



Figure 5. A plot of the flux F(x = 0) at spin location x — 0. The measured error 
bars are smaller than the point size (errors of order 10~ 7 ) and are not included in the 
graph. 
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Figure 6. We compare the exact result in with simulations for small values of the 
spin position x in (a). The coloured lines in (a) are included to guide the eye. In (b), 
we look at larger values of x and show the result of the simulations (black squares), 
the exact solution (red disks), and also the Bessel function approximation (coloured 
lines) we derived in (|73[) . The simulation results yielding flux values indistinguishable 
from zero are not included. 



We can also test our exact result in ( 152|) by numerically evaluating the integral 
in Mathematica. The results are shown in figure E^a). The simulations again agree 
perfectly with the exact result for the fluxes F(x) for all tested values of x and 7. 
Notice that for small values of 7 and large values of x, the magnitude of F(x) is quite 
small and indistinguishable from zero in the simulations. Thus, for these points we just 
show the numerically integrated exact result. Figure EJb) also includes approximation 
( I73p . Notice that the approximation agrees remarkably well with the simulation and 
exact result for small values of x and captures the cross-over between the exponential 
and power-law decay of F(x) as 7 — > 1. 

Finally, simulations can be used as a tool to better understand the behavior of 
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Figure 7. The plot markers are simulation results for ln|.F(x)| for three separate 
chain simulations with different pairs of values for jh,c- The different marker styles are 
associated with different sets of coloured line segments, with colours representing the 
7/j.c pairs shown in the plot legend. The solid line segments are best fits of the data in 
each bath for each simulated chain to the function In |F(ir)| = —\x\/£h,c ~ m M/2 + b, 
where b is our fitting parameter and is given by (|63|) with 7 = jh.c- For spins 
with x < 0, we increment a; by 1 in the fitting function order to avoid the logarithmic 
singularity at x = 0. The dashed magenta line simply connects the j c = 1 points as 
we do not have a conjectured fit function for this data. 

the model for arbitrary values of the bath temperatures and T c and to point us 
toward possible analytic solutions. For example, we can check to see if the asymptotics 
we found in f )62|) for T c > and — > 00 describe the flux decay when the hot bath 
is at a finite temperature. In figure El we fit simulation results for ln|F(x)| to the 
function —\x\/£(Th jC ) — (\n\x\)/2 + A, where £(7\ iC ) is half of the equilibrium correlation 
length for a chain at temperature c and A represents the fitting parameter and the 
overall amplitude of |_F(x)|. Thus, we find that the simulations are consistent with 
decaying as ~ Ae~^^ lyTh ^ / \/\x\ in both baths. This result confirms the 

intuitive picture of the flux in figure [U We see that the flux in the hot bath must decay 
faster than in the cold one since we will have ^(T^) < £(T C ). An analytic check of this 
conjecture would be a natural extension of the results presented here. In the T c — > 
(7 C — y 1) limit, the flux decays more slowly into the cold bath, as shown by the points 
connected by the dashed line in figure [71 However, we seem to lose the x~ 3 power law 
behavior for Th < 00 and find a different kind of long range decay. We have no conjecture 
for the behavior of this decay, so no fit to the simulation data was attempted. We also 
do not understand how the overall amplitude A depends on Th and T c . All we know is 
that it should decay to zero as — > T c and the whole chain approaches equilibrium. 
This is evident in figure [7J where the chains with values of 7^ close to 7 C have smaller 
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values of at the junction around x = 0. 

7. Conclusions and Outlook 

We have now examined in detail the exact results available for the energy flux through 
two Ising chains at different temperatures. This flux is an interesting quantity as it is 
a purely non-equilibrium quantity which vanishes at equilibrium. It is also physically 
relevant, as it describes the rate at which energy is transferred between the Ising chains 
due to a temperature gradient. We found that when one of the chains is at infinite 
temperature (T/, — > oo), an energy flux with magnitude u = 2/7 + a/4/7 2 — 1 with 
7 = tanh(2/3J) is injected into the cold chain from the hot chain. This energy then 
gets transported along the cold chain and eventually dissipates. We have given an 
exact expression for this flux F(x) and also an excellent approximation. The main 
features of this flux is an exponential decay for T c > and large x with a decay length 
£ = —{In [tanh(/3 J)]} _1 /2 equal to half of the correlation length of an infinite chain at 
equilibrium at temperature T c . As T c — > 0, we get the maximum amount of energy 
injected into the cold chain and the flux F(x) decays much more slowly, as the power 
law x~ 3 at T c = 0. 

Another way of understanding the T c — > limit is by focusing on the dynamics of 
the boundaries between clusters of up and down spins. We can treat these boundaries 
as particles by identifying any two spins that are anti-aligned. For example, we can 
have spin configurations tt 44 or 44 T1\ where we use the o symbol to denote the 
boundary particle. Notice that it costs no energy to flip either the second or third spin. 
Thus, when either of these spins is selected in a Monte Carlo simulation, it will flip with 
a probability of 1/2 , moving the boundary either to the left or to the right. Thus, these 
domain boundaries will perform unbiased random walks along the chain coupled to the 
cold bath with T c = until they are removed at the two ends of the chain or until two 
domain boundaries annihilate via the transitions | I T — >TTT an d 4 t 4~ ^444- 
This pair-annihilation will release an energy of 4J into the cold bath. Notice that the 
infinite temperature x = spin will generate these particles at one end of the chain as 
it is allowed to flip even when its aligned with its neighbor at x — 1. Since this is the 
only mechanism by which we change the energy of the cold spin chain, when 7 = 1, 
F(x) measures the frequency of these annihilation events along the chain due to the 
injection of these boundary "particles" at x — 0. These kinds of particle models (and 
their Ising model duals) have been the subject of much study p2J [13], ED, E2] and many 
exact results for the steady state distributions and correlation functions are available via 
fermionic methods, Bethe ansatz techniques, etc. This dual language might be useful 
in understanding what happens in the cold bath when T c = and Th < 00, as we were 
only able to fully understand what happens in the limit T h — > 00. When T c > 0, F(x) 
also includes contributions from pair creation of the boundary particles. 

A natural generalization of the model considered in this paper to two dimensions 
is to look at spins on a square lattice. The lattice can be partitioned into two halves, 
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with spins in the two halves coupled to two different heat baths with temperatures Th 
and T c . In two dimensions, there is a phase transition at finite temperature and the 
ordered phase (aligned spins) becomes stable at low temperatures. This means that at 
low temperature, more energy is required to both form and move domain boundaries. 
The boundaries at low temperature do not move as freely inside the system as they 
do in one dimension because spins at the boundaries will be typically surrounded by 
more aligned than anti-aligned spins. Thus, we expect the fluctuations at the interface 
between the two halves of the lattice not to propagate very far at low temperatures. 
The total energy flux between hot and cold baths coupled to spins on a two-dimensional 
lattice was studied previously [231 El] to test the validity of fluctuation relations out 
of equilibrium. It would be interesting to expand upon this research by looking at the 
spatial dependence of the flux. For example, one could study the generalization of F(x): 
the average rate of energy change due to spin flipping at a distance x away from the 
interface. 

Preliminary simulation results on a square lattice (with periodic BCs along the 
direction parallel to the interface between the two lattice halves and open BCs on the 
other two sides) suggest that, unlike in one dimension, the flux F(x) is vanishingly small 
when T/j = oo and T c = 0. The flux decays approximately exponentially for T h = oo 
and arbitrary finite T c > 0. We find that, assuming an exponential decay of F(x), 
the associated decay length is small for all T c (on the order of the lattice spacing) and 
increases as T c approaches the critical temperature T* ~ 2.3 J [18] for both T c < T* and 
T c > T*. Also unlike the one- dimensional case, the decay length of F(x) does not appear 
to be simply related to the equilibrium correlation length. Of course, these results need 
to be checked with more careful simulations, especially when T c is close to T*, where 
there might be novel critical behavior and associated power laws. 

Our simulations were limited to small system sizes (squares of about 20 spins on a 
side) and we expect that there are finite-size effects. A comprehensive account of these 
effects is beyond the scope of this paper, and the two-dimensional generalization of the 
model remains an interesting open problem for future study. Such a study would require 
a detailed look at all the possible combinations of Th and T c relative to T*, as each half 
of the lattice is expected to behave qualitatively differently depending on whether its 
temperature and that of its partner is above or below the critical temperature. As 
mentioned in the previous paper [T3], we also expect a nontrivial magnetization profile 
at temperatures below the critical temperature. Finally, the preliminary simulation 
results suggest that an interesting case is T h = oo and T c = T*, where we seem to have 
the largest energy flux. 

There are some more obvious open problems, like the analytic behavior of F(x) in 
one dimension for the case where both T c and Th are finite. We have formal expressions 
for this case, but we do not understand precisely how F(x) decays away from the junction 
for this more general case. We conjecture in section [6] that for non-zero temperatures 
F(x) decays exponentially into both baths. This should be checked rigorously. As 
mentioned in section [6], we also do not know the precise way in which F(x) — > 
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as the entire chain moves closer to equilibrium as T c — > T^. Finally, as discussed in 
[H], there are countless ways to extend the calculations we presented here. We could 
include different dynamics (such as Kawasaki spin-exchange), explore different boundary 
conditions, etc. 

It is interesting that the far-from-equilibrium system studied here displays behavior 
reminiscent of a second order phase transition when — > oo and T c — > 0: We have a 
diverging decay length and a transition to a power law. However, unlike an equilibrium 
phase transition, we do not have a general framework for understanding this crossover. 
This means that exact results such as the ones presented in this work are especially 
useful in guiding more general studies of non-equilibrium phenomena. 
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